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ABSTRACT 


A  method  for  attacking  a  wide  variety  of  commonly  occurring  heat- 
transfer  problems  is  presented.  The  method  stems  from  a  mathematical 
model  obtained  by  dividing  a  system  into  pieces.  The  differential  equa¬ 
tion  of  the  temperature  of  each  piece  with  respect  to  time  is  derived 
considering  heat  sources  and  heat  transfer  by  conduction  and  radiation. 
The  resulting  system  of  first-order  ordinary  differential  equations  is 
put  in  a  form  convenient  for  programming  a  numerical  solution  by 
defining  coefficients  in  terms  of  the  data  of  the  system.  Form-surface 
factors  are  a  generalization  of  form  factors  to  include  multiple  reflec¬ 
tions.  The  factor  is  defined  and  discussed.  A  formula  is  derived  for 
the  factor  in  terms  of  the  form  factors  and  the  reflectivities.  The  form- 
surface  factor  relations  for  the  conservation  of  energy  and  the  reciproc¬ 
ity  law  are  derived.  The  conditions  necessary  and  sufficient  for  its 
existence  are  derived.  Methods  of  calculation  are  discussed.  A 
FORTRAN  subroutine  is  presented  that  calculates  the  form-surface 
factors  from  the  form  factors  and  the  reflectivities.  A  FORTRAN  sub¬ 
routine  is  presented  which  performs  unit  conversions  on  data  tabulated 
for  the  mathematical  model.  Another  FORTRAN  subroutine  is  presented 
which  calculates  the  coefficients  from  the  tabulated  data.  A  numerical 
method  for  solving  the  system  of  differential  equations  is  developed  and 
a  FORTRAN  subroutine  is  presented  which  effects  the  solution.  Several 
example  problems  are  worked  out,  and  the  results  are  compared  with 
analytical  solutions. 
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SECTION  I 
INTRODUCTION 


General  methods  are  needed  to  solve  wide  varieties  of  commonly 
occurring  heat-transfer  problems  to  shorten  the  time  from  defining  the 
problem  to  obtaining  a  solution..  This  report  presents  a  method  for 
transient,  conduction,  and  radiation  heat -transfer  problems.  The 
given  system  is  divided  into  pieces  to  form  a  mathematical  model. 

The  differential  equation  of  each  piece  with  respect  to  time  is  derived 
considering  internal  power  sources  and  heat  transfer  by  conduction  and 
radiation.  The  resulting  system  of  first-order  differential  equations  is 
put  in  a  form  convenient  for  programming  a  numerical  solution  by  de¬ 
fining  coefficients  in  terms  of  the  data  of  the  system.  FORTRAN  sub¬ 
routines  are  presented  which  convert  units  of  the  data  of  the  system, 
calculate  the  coefficients,  and  effect  a  numerical  solution  to  the  equa¬ 
tions. 

The  form-surface  factor  is  a  generalization  of  the  form  factor  to 
include  multiple  reflections.  This  factor  is  defined  and  discussed,  and 
a  formula  is  derived  giving  the  factor  in  terms  of  the  form  factors  and 
the  reflectivities  of  the  system.  The  form-surface  factor  relations  for 
conservation  of  energy  and  the  reciprocity  law  are  derived.  Also,  con¬ 
ditions  necessary  and  sufficient  for  the  existence  of  the  form-surface 
factor  are  derived.  Methods  of  calculation  are  discussed  and  a 
FORTRAN  subroutine  is  presented  which  calculates  the  form-surface 
factor  from  the  form  factors  and  the  reflectivities. 

Once  familiar  with  the  mathematical  model  and  the  subroutines,  the 
following  steps  are  taken:  (1)  The  system  is  divided  into  pieces;  (2)  a 
FORTRAN  program  is  written  to  calculate  and  tabulate  data  for  the 
mathematical  model;  (3)  a  subroutine  is  called  which  calculates  the 
form-surface  factors  from  the  form  factors  and  the  reflectivities;  (4)  a 
subroutine  is  called  which  performs  unit  conversions  on  all  the  tabu¬ 
lated  data;  (5)  a  subroutine  is  called  which  calculates  the  coefficients 
from  the  converted  data;  (6)  lastly,  a  subroutine  is  called  which  effects 
the  numerical  solution  to  the  system  of  differential  equations. 

Example  problems  illustrating  the  method  are  worked  out,  and  the 
results  are  presented  and  compared  with  analytical  solutions.  The 
FORTRAN  listings  of  the  programs  used  in  the  example  problems  are 
given  in  appendixes  to  be  used  as  a  guide. 
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SECTION  II 

MATHEMATICAL  MODEL 


2.1  DERIVATION  OF  EQUATIONS 

The  data  input  to  the  computer  programs  are  based  on  a  mathe¬ 
matical  model  of  the  system.  The  model  is  obtained  by  dividing  the 
system  into  pieces,  with  uniform  thermal  properties  (Fig.  1).  The 
equations  for  this  mathematical  model  are  derived  by  performing  an 
energy  balance  on  each  piece. 


Fig.  1  Illustrative  System 
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Assume  that  the  model  consists  of  n  pieces.  The  rate  of  tempera¬ 
ture  change  of  the  ith  piece  depends  on  the  net  rate  that  heat  is  gained 
by  the  piece  according  to  the  equation 

di  Vi  cs  T i  =  Pj  ( 1) 

where  the  dot  represents  the  time  derivative.  The  types  of  heat  trans¬ 
fer  which  will  be  considered  are  heat  sources,  conduction,  and^radiant 
heat  transfer  in  a  vacuum;  thus,  the  net  rate  that  heat  is  gained  is 

Pi  =  Pj-  -  Psfc  +  (Pi-  -  Pi*)  (2) 

The  first  term  accounts  for  heat  supplied  by  a  source  such  as  an  elec¬ 
trical  heater.  This  quantity  may  be  negative  in  which  case  it  represents 
a  heat  sink. 

The  rate  that  heat  is  gained  by  the  ith  piece  by  conduction  is 

P>=  Ik/  AT  •  da]  (3) 

joniL  Mj  J 

where  the  summation  is  over  the  pieces  j  in  contact  with  the  ith  piece 
(abbreviated  ]  on  i)  and  where  the  differential  area  vector  is  directed 
outward.  It  will  be  assumed  that  the  temperature  gradient  is  constant 
over  each  conducting  area,  ajj,  and  that  it  is  proportional  to  the  tem¬ 
perature  difference  Tj  -  Ti.  To  derive  the  proportionality  constant 
consider  the  linearized  case  illustrated  in  Fig.  2. 


T 
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It  can  be  seen  that  there  will  be  a  discontinuity  of  the  temperature 
gradient  at  the  interface  of  the  two  pieces  if  their  thermal  conductivities 
are  different.  Conservation  of  energy  requires  that  at  the  interface 


(4) 


Using  this  relation  and  geometrical  reasoning  one  derives 

dT|  2(Tj-  -  Tf) 


(5) 


Using  these  approximations  Eq.  (3)  becomes 

pi  _  2  2  *i  *j  aij  (Tj  -  Tj) 

i  on  i  djj  (#<£  +  Kj)  ^ 

where  the  summation  is  over  the  pieces  which  are  in  contact  with  the 
ith  piece.  Note  that  the  conducting  area,  ay,  is  defined  as  the  area 

perpendicular  to  the  conducting  path,  thus  the  cosine  factor  resulting 
from  the  dot  product  of  Eq.  (3)  is  one.  The  conducting  distance,  dy, 
is  the  distance  along  the  conducting  path. 


The  rate  that  radiant  energy  is  emitted  by  the  ith  piece  is 

Pi"  =  S;  fi  a  Ti4  (7) 


The  rate  that  radiant  energy  is  absorbed  by  the  ith  piece  is 

Pi8  =  ji,  GjiSj-ejvTj" 


(8) 


The  factor  Gji  is  defined  as  the  ratio  of  that  part  of  the  flux  incident  on 

the  ith  piece,  after  all  reflections,  which  was  originally  emitted  by  the 
jth  piece  to  the  flux  emitted  by  the  jth  piece.  Thus  the  last  four  factors 
of  Eq.  (8),  yielding  the  rate  at  which  radiant  energy  is  emitted  by  the 
jth  piece,  are  multiplied  by  Gji  to  obtain  the  resulting  flux  on  the  ith 
piece.  This  is  multiplied  by  the  absorptivity  and  summed  over  all  the 
pieces  to  get  the  rate  heat  is  gained  by  absorption  by  the  ith  piece.  The 
explicit  dependence  of  Gji  on  the  geometry  and  surface  properties  of  the 
system  is  derived  in  the  next  section  along  with  some  important  relations 
it  must  satisfy.  A  FORTRAN  subroutine  is  described  in  Section  IV  which 
can  be  used  for  its  calculation. 
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After  all  substitutions,  the  energy  balance  equation  for  the  ith  piece, 
Eq.  (1),  becomes 


*l  Vi  ci  Ti  =  Pi8 


2Kj  Kj  ajj  (Tj  —  T{) 
+  i  ^ 

j  o  n  i  d I j  (Kj  4  Kj) 


Si  fi^Ti4  +  21  ai  Gji 


Sj  fj  a  Tj4 


(9) 


There  is  one  equation  like  Eq.  (9)  for  each  piece  in  the  model. 


Equation  (9)  can  be  written  in  a  more  convenient  form  for  pro¬ 
gramming  a  numerical  solution  by  collecting  terms  and  solving  for  Ti: 


where 


Ti  =  Ai  +  2BijTj+  ^CijTj4 


Ai 


P" 

dj  Vi 


Bij 


0 


2  «1  *j  «ij _ 

V{C|  dij  (Kj  4  Kj) 

2*i  j.  KkBik 
di  Vi  Ci  koni  djfc  (Ki  4  Kfc) 


j  not  on  i 
j  on  i 


j  =  > 


CiJ 


a i  G]l  Sj  fj  a 

di  Vi  ci 


Sj  fj  g 

di  W  d 


(1  -  Q.j  Gjj) 


j 


‘  =  j 


(10) 

(11) 


(12) 


(13) 


A  FORTRAN  subroutine  is  described  in  Section  V  which,  when  given 
the  values  of  the  quantities  on  the  right  side  of  Eqs.  (11)  through  (13), 
will  evaluate  these  equations.  When  evaluating  Eq.  (12),  the  subroutine 
assumes  that  pieces  are  in  contact  if  and  only  if  their  conducting  area 
is  positive.  The  results  of  this  subroutine  can  be  used  directly  in 
another  subroutine  described  in  Section  VI  which  numerically  calcu¬ 
lates  the  solution  of  the  system  of  differential  equations  represented 
by  Eq.  (10). 


Dividing  the  system  into  hypothetical  pieces  has  effectively  reduced 
the  partial  differential  integral  equation  relating  temperature  as  a  func¬ 
tion  of  time  and  position  into  a  system  of  ordinary  differential  equations 
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for  which  there  are  standard  general  numerical  methods  (Ref.  1)  for 
solution.  This  method  of  reducing  partial  differential  equations  to 
ordinary  differential  equations  is  discussed  on  page  110  of  Ref.  1  which 
lists  additional  references. 


2.2  BOUNDARY  CONDITIONS 

Most  heat-transfer  problems  are  subject  to  boundary  conditions, 
usually  in  the  form  of  an  insulated  boundary  or  a  boundary  held  at  a 
fixed  temperature.  The  insulated  boundary  requires  that  no  heat  flows 
across  the  boundary.  Since  the  mathematical  model  was  based  on  an 
isolated  system,  this  type  of  condition  is  inherent  to  the  model  and  re¬ 
quires  no  further  analysis.  The  fixed  temperature  boundary  condition 
can  be  imposed  by  including  the  boundary  as  one  of  the  pieces  of  the 
model,  say  the  ith  piece,  and  letting  its  total  heat  capacity,  d^  cj, 
be  infinite.  Thus,  no  amount  of  heat  transferred  into  or  out  of  the 
piece  will  change  its  temperature.  This  may  or  may  not  be  a  true 
representation  of  the  actual  physical  system,  but  it  gives  the  desired 
result  of  constant  temperature.  Mathematically,  from  Eqs.  (11) 
through  (13),  it  is  seen  that  when  the  total  heat  capacity  of  the  ith  piece 
is  infinite,  then 

A;  =  0 

Bjj  -  0  j  -  1,  2,  .  .  .  ,  n 

C  j j  =  0  j  =  1,  2,  .  .  .  ,  n 

thus  Eq.  (1)  becomes 

T;  =  0 

so  the  temperature  is  constant.  The  FORTRAN  subroutine  described 
in  Section  V  was  programmed  to  recognize  this  case  and  set  the  coeffi¬ 
cients  of  the  ith  equation  equal  to  zero  when  the  input,  ci,  is  zero. 


SECTION  III 

FORM-SURFACE  FACTORS 


3.1  DEFINITION 

Let  Qij  be  the  portion  of  the  flux  incident  on  the  jth  piece  which  was 
originally  emitted  by  the  ith  piece.  Let  Wj  be  the  flux  emitted  by  the 
ith  piece,  for  radiant  flux  that  is 

W;  =  S;  q  <7  T j4  (14) 
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The  form-surface  factor  is  defined  as 


Qij 

Wi 


(15) 


This  definition  warrants  some  elaboration.  It  is  tacitly  assumed  that 
flux  can  be  distinguished  as  to  its  origin.  Picture,  for  instance, 
radiant  energy  as  a  stream  of  photons.  Let  the  units  of  flux  be  photons 
per  second  and  assume  that  the  average  time  from  a  photo  being  emitted 
to  it  being  absorbed  is  negligible.  Consider  a  photon  which  is  emitted 
by  the  ith  piece,  is  reflected  around  the  system  striking  the  jth  piece 
six  times,  and  finally  is  absorbed  by  the  kth  piece.  This  photon  will 
contribute  one  unit  to  Wi  and  six  units  to  Q^j.  Statistically,  as  the  num¬ 
ber  of  photons  becomes  large,  the  ratio,  Eq.  (15),  will  become  constant. 
It  can  be  seen  that  depends  only  on  the  geometry  and  surface  prop¬ 
erties  of  the  system;  thus  its  name:  the  form-surface  factor.  With  this 
idea  in  mind  it  is  no  longer  necessary  that  flux  be  distinguished  as  to 
its  origin.  It  is  only  necessary  to  be  able  to  associate  with  the  flux 
emitted  by  the  ith  piece  a  quantity  of  flux,  Qy,  incident  on  the  jth  piece, 

for  which  the  ratio,  Eq.  (15),  will  be  constant  as  W*  is  varied.  Ob¬ 
viously,  the  ratio  of  Eq.  (15)  will  depend  on  the  distribution  of  Wi  over 
the  ith  piece  and  also  on  the  angular  emittance  at  each  point.  The  form- 
surface  factor  will  be  defined  by  Eq.  (15)  where  Wi  is  emitted  diffusely 
and  uniformly. 


3.2  DERIVATIONS 


There  are  several  ways  to  derive  an  equation  for  the  matrix,  G,  of 
the  form-surface  factors  of  a  system.  Three  derivations  are  included 
here  because  they  lead  to  a  better  understanding  of  the  factor  and  be¬ 
cause  some  interesting  observations  can  be  made  on  the  respective 
approximations . 


Consider  the  integral  equation  for  the  irradiation 


H«)  =  /p(?>H(?)F(?,jf)dSy  -  Jf(?)<rT4<?)  F(?,  Z)  dSy  ( 16) 

s  s 

where 


-*  ,  FdSy-ASx 

F(y,  x)  =  lim  - — -  (17) 

ASx-»0 

where  ASX  is  a  small  area  around  the  point  x.  If  the  integrals  are 
approximated  using  the  partitioning  as  determined  by  the  mathematical 
model  as  described  in  Section  II,  one  obtains 


H: 


n 

=  2 


i=  1 


p)Hj 


Fji 


Si 


n 

2 


(,  aT,- 


Si 


(18) 
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Multiplying  through  by  Si  this  becomes 

Qi  =  PjQj  F,i  +  |iWj  Fji  (19) 


Letting  Q  and  W  be  column  matrices  of  the  Qi  and  Wi,  respectively, 
and  letting  R  be  the  diagonal  matrix  of  the  reflectivities,  Eq.  (19)  can 
be  written  as  the  matrix  equation 


Q  =  F'RQ  +  F'W 

(20) 

Solving  for  Q  one  obtains 

Q  =  (I  -  F 'R)~‘  F'W 

(21) 

Defining 

G  =  F(I  -  RF)-* 

(22) 

Equation  (21)  becomes 

Q  =  G'W 

(23) 

or 

n 

Qj  =  i£i  Gii  Wi 

(24) 

From  Eq.  (24)  it  is  seen  that  a  quantity  of  flux,  Q^,  incident  on  the 

jth  piece  given  by 

Qij  =■  G;j  Wj 

(25) 

can  be  associated  with  the  flux  emitted  by  the  ith  piece. 

From  Eq.  (22) 

it  is  seen  that  G  depends  on  F  and  R  which  are  determined  by  the  geom¬ 
etry  of  the  system  and  its  surface  properties,  respectively.  Further¬ 
more,  G  does  not  depend  on  W;  thus  the  elements  Gij  of  the  matrix,  G, 
conform  to  the  definition  of  the  form-surface  factor. 

Equation  (22)  can  also  be  derived  by  summing  all  the  reflections. 

Let  Qk  be  the  column  matrix  of  the  incident  flux  which  has  been  reflected 
exactly  k  times.  The  flux  incident  directly  is 

Q°  =  F'W  (26) 

If  it  is  assumed  that  flux  is  redistributed  uniformly  over  a  piece  before 
being  reflected,  then  the  flux  vector  for  the  next  reflection  can  be  ob¬ 
tained  from  the  previous  one  by  the  equation 

Qk  =  f'rQ1"1  k  >  0  (27) 
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The  total  incident  flux  is  the  sum  of  all  the  reflections  or 

Q  =  S  Qk 

k=0 

=  F'W  +  F'RQ0  +  F'RQ1  +  .  .  . 

=  F'W  -  F'R(Q°  +  Ql  -  ...) 
or 

Q  =  F'W  +  F'RQ 


(28) 

(29) 


This  is  the  same  as  Eq.  (20);  thus  Eq.  (22)  for  G  follows  the  same  as 
before. 


A  different  equation  for  G  can  be  obtained  from  an  equation  that  is 
consequent  from  the  definition  of  G.  The  flux  incident  on  the  jth  piece 
resulting  from  is  the  sum  of  the  flux  incident  directly  from  the  ith 
piece  and  the  flux  from  all  reflections.  If  it  is  assumed  that  flux  inci¬ 
dent  on  a  piece  directly  from  the  ith  piece  is  redistributed  uniformly 
over  the  piece  before  being  reflected  the  first  time,  then  the  form- 
surface  factor  can  be  used  to  obtain  the  flux  incident  on  the  jth  piece 
from  all  reflections.  Thus  it  is  deduced  that 

Qij  =  Wj  Fij  +  k^l  Wi  F;k  Pic  Gkj  (3Q) 

Thus,  from  Eq.  (15)  is  obtained 

n 

Gij  “  Fij  +k|lFikPkGkj  (3!) 


from  which  is  obtained  the  matrix  equation 

G  =  F  -  FRG  (32) 

Solving  for  G  one  obtains 

G  =  (I  -  FR)"1  F  (33) 

Equations  (22)  and  (33)  are  two  different  equations;  however,  they  give 
the  same  matrix,  G,  as  is  seen  from 

F(I  -  RF)-  =  (I  -  FR)"‘F  <34) 

Multiplying  each  side  of  this  equation  on  the  left  by  (I  -  FR)  and  on  the 
right  by  (I  -  RF)  one  obtains 

* 

(I  -  FR)  F  (I  -  RF)"1  (I  -  RF)  =  (I  -  FR)  (I  -  FR)-‘F  (I  -  RF) 
or 

(I-FR)F  =  F(I-RF)  (35) 
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This  multiplication  does  not  change  the  equality  since  the  multipliers 
are  nonsingular.  Equation  (35)  is  indeed  an  equality  since  by  the 
distributive  law  it  is  equivalent  to 

F  -  FRF  =  F  -  FRF  (36) 

which  verifies  the  equivalence  of  Eqs.  (22)  and  (33). 

*  It  is  interesting  to  note  that  an  equation  for  G  has  been  derived  in 
three  different  ways  using  three  different  simplifying  approximations. 
The  approximation  in  the  first  derivation  was  the  discretization  of  the 
integrals.  In  the  second  derivation  flux  was  redistributed  uniformly 
over  a  piece  before  each  reflection,  whereas  in  the  third  derivation  the 
flux  was  redistributed  uniformly  over  each  piece  before  only  the  first 
reflection.  Nevertheless,  the  matrices  obtained  for  G  in  all  three 
derivations  are  identical. 


3.3  CONSERVATION  OF  ENERGY 

From  the  definition  of  the  form  factor  one  obtains  for  a  closed  sys¬ 
tem  the  relation  that  the  sum  of  each  row  of  F  must  be  one,  or 

JjFij  =  1  (37) 

Since  the  form-surface  factor  is  a  generalization  of  the  form  factor, 
one  would  look  for  a  relation  for  G  generalizing  Eq.  (37).  This 
generalization  is  easily  obtained  from  the  definition  of  the  form-surface 
factor  and  the  conservation  of  energy.  Since  in  a  closed  system,  the 
energy  emitted  by  a  piece  must  all  be  absorbed, 

n 

E  Wj  Gjj  at  =  Wj 
i=l  J 

or  dividing  by  Wj 

D 

E  a\  Gjj  =  1 

j=l 

Thus,  from  physical  laws,  the  form-surface  factors  must  satisfy 
Eq.  (39).  However,  it  does  not  necessarily  follow  that  the  matrix,  G, 
given  by  Eq.  (22)  satisfies  Eq.  (39),  since  Eq.  (22)  was  derived  under 
approximate  conditions.  If  Eq,  (22)  is  to  be  used  to  calculate  G,  it  is 
imperative  that  it  be  known  how  well  the  resulting  G  satisfies  Eq.  (39), 
since  conservation  of  energy  in  the  mathematical  model  depends  on 
Eq.  (39).  Fortunately,  it  can  be  proved  mathematically  that  the  G 
given  by  Eq.  (22)  satisfies  Eq.  (39)  exactly. 


(38) 

(39) 
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Multiplying  both  sides  of  Eq.  (22)  on  the  right  by  (I  -  RF)  one  ob¬ 
tains 

G  -  GRF  =  F  (40) 


Summing  row  i  of  Eq.  (40)  one  obtains 


£  Gij  -  £  £  G it  PkFkj 

j=i  J  i=  1 i 


(41) 


Rearranging  this  becomes 


(42) 


But,  by  Eq.  (37),  the  quantities  in  parentheses  are  one,  so  Eq.  (42) 
becomes 

n  n 

.*Gji  ~kl,Gikpk  =  1  (43) 


Since  k  in  the  second  summation  is  a  dummy  index  it  can  be  changed 
to  j  and  the  two  summations  can  be  combined.  Factorization  of  Gy 
and  substitution  of  orj  for  (1  -  Pj)  gives  the  desired  result,  Eq.  (39). 


3.4  RECIPROCITY  LAW 

It  can  be  proved  that  form  factors  obey  the  reciprocity  law,  that  is 
that  SF  is  symmetric  or 

SF  =  (SF)'  (44) 

where  S  is  the  diagonal  matrix  of  the  S^.  It  can  be  proved  from  the 
second  law  of  thermodynamics  that  the  form -surface  factor  also  obeys 
the  reciprocity  law.  From  Eq.  (8)  the  net  rate  of  heat  transfer  from 
the  ith  piece  to  the  jth  piece  is 

Pjj  =  aj  Gij  Si  q  o  Tj4  -  a,  Gjj  Sj  fj  a  Tj4  (45) 

If  Ti  =  Tj  =  To  then  the  second  law  of  thermodynamics  requires  that 

Q  “ 

Py  be  zero.  Setting  the  emissivities  equal  to  the  absorptivities  this 
condition  becomes 

dj  dj  a  Tq4  (Si  Gij  -  Sj  Gji)  =  0  (46) 

In  general  the  factors  on  the  left  are  not  zero,  thus  the  quantities  in 
parentheses  must  be  zero.  Thus  SG  must  be  symmetric  or 

SG  =  (SG)'  (47) 
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By  reasoning  parallel  to  that  in  Section  3.  3,  it  is  desirable  to  prove 
mathematically  that  the  matrix,  G,  given  by  Eq.  (22),  satisfies  Eq.  (47); 
since  this  is  required  if  the  mathematical  model  is  to  obey  the  second 
law  of  thermodynamics . 

Multiplying  both  sides  of  Eq.  (22)  on  the  left  by  S  one  obtains 

SG  =  SF  (I  -  RF)-1  (48) 

Using  the  substitution,  Eq.  (44)  for  SF,  and  the  fact  that  S  is  diagonal, 
thus  symmetric,  Eq.  (48)  becomes 

SG  =  F'S(I  -  RF)-1  (49) 

Since  S  is  diagonal  and  has  no  zero  elements  on  the  diagonal,  it  is  non¬ 
singular  and  one  can  write 

SG  =  F '  [(I  -  RF)  S-1]-1 
=  F'[S-1  -  RFS-‘]-‘ 

=  F  '[S-1  -  S-1  SRFS-1]-*  (50) 

Since  S  and  R  are  both  diagonal,  they  commute;  thus 

SG  =  FIS"1  -  S“* RtSFjS-*]--1  (51) 

Making  the  substitution  of  Eq.  (44)  for  SF  again  one  obtains 

SG  =  F'[S“*  -  S^RF'SS-*]-* 

=  F'[S~l  -  S"‘  RF']"* 

=  F  '[S-1  (I  -  RF')]-1 
=  F'[I  -  RF']-1  S 

=  [S  (I  -  FR)-1  FJ'  (52) 

Thus  by  Eq.  (33),  which  has  been  proven  to  yield  the  matrix  G  identical 
to  Eq.  (22),  Eq.  (52)  becomes 

SG  =  (SG)' 

and  it  has  been  proved  that  the  G  given  by  Eq.  (22)  obeys  the  reciprocity 
law. 


3.5  SINGULARITY  AND  NONSINGULARITY  CONDITIONS 

Equation  ( 22)  is  a  formula  for  G  provided  the  matrix  (I  -  RF)  is  non¬ 
singular.  It  is  desirable  to  know  the  general  conditions  for  the  singularity 
or  nonsingularity  of  this  matrix.  These  conditions  can  be  induced  from 
physical  considerations  and  then  proved  mathematically.  Consider  each 
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subsystem  of  a  given  system  which  is  by  itself  a  closed  system.  If  for 
a  given  subsystem  all  the  reflectivities  are  one,  then  any  flux  inside  the 
subsystem  will  be  reflected  continuously  and  G  for  the  system  will  not 
exist.  If,  on  the  other  hand,  every  subsystem  has  at  least  one  surface 
with  a  reflectivity  less  than  one,  then  any  flux  in  the  system,  eventually 
after  multiple  reflections,  will  be  absorbed  and  G  will  exist.  These 
conditions  can  be  stated  as  follows:  Group  the  pieces  into  sets 
Il»  l2»  •  •  •  Im  such  that  for  each  set  Ifc 


Then  by  renumbering  the  pieces  the  form  factor  matrix  can  be  parti¬ 
tioned  into  a  diagonal  matrix,  the  elements  of  which  are  the  subform 
factor  matrices  of  the  groups  described  above.  The  matrix  (I  -  RF) 
for  the  system  is  nonsingular  if  and  only  if  the  matrix  (I  -  RF)  for  each 
subsystem  is  nonsingular.  From  the  above  physical  considerations,  a 
matrix  for  a  subsystem  is  nonsingular  if  and  only  if  at  least  one  of  its 
reflectivities  is  less  than  one.  The  property  of  the  subform  factor 
matrices  described  above  is  seen  to  correspond  to  the  mathematical 
term,  irreducible.  A  general  matrix  F  is  said  to  be  reducible  if  two 
sets  of  subscripts  I  and  J  exist  such  that: 

1.  There  exists  i  such  that  i  e  I 

2.  There  exists  j  such  that  j  €  J 

3.  k  e  I  U  J  for  every  k 

4.  k  4  I  J1  J  for  every  k 

5.  If  i  e  I  and  j  e  J,  then  Fy  =  0 

Otherwise  the  matrix  is  irreducible. 


If  (I  -  RF)  is  singular  then  there  exists  a  nontrivial  solution  to 

(I  -  RF)  X  =  0  (53) 

In  terms  of  the  elements  this  becomes 

xi  -  Pi  £  Fijxj  =  0 
j=i  J 

or 

xi  “  Pi  Fii  xi  =  Pi  J.  Fij  xj 

J=»> 


Let  x'  be  the  value  of  the  largest  element  in  X,  then 


or 


xi  —  Pi  Fjj  xj  <  pi  x  X  Fji 

j=/i 

xi  -  Pi  F ii  xi  <  pi  x'U  -  Fh) 


(54) 
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In  particular  if  =  x',  then  Eq.  (54)  becomes  for  i  =  k 

xic  -  p k  Fkk  *k  £  Pk  xk  -  Pk  Fkk  *k  (55) 

The  only  possible  way  for  Eq.  (55)  to  hold  is  for  Pk  to  equal  one  and  for 
the  equality  of  Eq.  (54)  to  hold.  Thus  for  (I  -  RF)  to  be  singular  the 
equality  of  Eq.  (54)  must  hold  for  i  equal  to  k.  That  is 

xj  =  x'  for  Fkj  4  0 

Then,  however,  Eq.  (55)  must  hold  for  every  k  such  that 

xk  =  x' 

The  only  possibility  that  there  could  be  elements  of  X  not  equal  to  x'  is 
for  there  to  exist  two  sets  of  subscripts  I  and  J  such  that  if  k  €  I  and 
j  e  J  then  Fkj  =  0,  that  is  for  F  to  be  reducible.  If  k  e  I,  then  Xk  must 
equal  x';  thus  Eq.  (55)  must  hold,  which  requires  pk  to  be  one.  Thus  a 
necessary  and  sufficient  condition  for  (I  -  RF)  to  be  singular  is  for  the 
system  itself  to  be  irreducible  and  all  of  its  reflectivities  be  one,  or 
for  there  to  exist  an  irreducible  subsystem  for  which  all  the  reflectivi¬ 
ties  are  one.  If  every  irreducible  subsystem  has  at  least  one  surface 
with  a  reflectivity  less  than  one,  then  (I  -  RF)  will  be  nonsingular. 


SECTION  IV 

CALCULATION  OF  FORM-SURFACE  FACTORS 


Just  as  there  were  a  variety  of  ways  to  derive  the  formula  for  G, 
there  are  a  variety  of  ways  to  calculate  it:  namely,  matrix  inversion, 
summation,  and  iteration.  The  obvious  way  is  directly  from  the 
formula,  Eq.  (22),  using  a  matrix  inversion  routine.  This  has  the 
advantage  of  being  simple,  and  since  it  is  direct,  no  convergence 
criterion  is  needed.  Though  it  was  not  apparent  at  the  outset,  matrix 
inversion  is  also  the  quickest  method.  However,  the  other  methods  will 
also  be  given  here,  since  it  is  sometimes  desirable  to  do  a  reflection 
by  reflection  analysis  of  a  system,  in  which  case  the  summation  and 
iteration  techniques  are  applicable,  whereas  the  matrix  inversion 
method  is  not. 


4.1  MATRIX  INVERSION  METHOD 

The  subroutine  GINV  was  written  in  FORTRAN  IV  to  calculate  G 
using  the  matrix  inversion  method.  The  routine  is  called  by  the  state¬ 
ment 

CALL  GINV  (N,  NDLM,  R,  F,  G,  D£T) 
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where 

N  Contains  n 

•  NDIM  Contains  the  dimension  of  F  and  G 

(see  example  below) 

R  Is  a  one -dimensional  array  containing 

the  reflectivities 

F  Is  a  two-dimensional  array  containing 

the  form  factors 

G  Is  a  two-dimensional  array,  in  which 

the  form-surface  factors  are  returned 

DET  Returns  the  value  of  the  determinant 

of  (I-RF)  which  will  be  zero  if  the 
routine  was  unable  to  invert  the 
matrix. 

As  an  example  of  the  purpose  of  NDIM,  suppose  the  main  program 
has  the  dimension  statement 

DIMENSION  R(60),  F(60,60),  G(60,60) 

If  G  is  to  be  calculated  for  a  system  of  ten  pieces,  then  the  call  state¬ 
ment  would  be 


CALL  GINV  (10,60,R,F,G,DET) 

Thus,  it  is  seen,  the  dimension  of  the  FORTRAN  variables  may  be 
greater  than  or  equal  to  (not  less  than)  the  order  of  the  algebraic 
matrices 

The  listing  of  GINV  is  given  in  Appendix  I.  The  subroutine  GINV 
calls  another  subroutine  INVERT  which  is  a  matrix  inversion  routine. 
Appendix  II.  The  timing  for  GINV  is  approximately 

t  =  t0n2  (3n  +  8) 

where  for  the  IBM  System/360  Model  50 

t0  =  1.43  x  10-4  sec 

It  is  to  be  mentioned  that  the  subroutines  set  aside  400  bytes  of  scratch 
storage  in  labeled  COMMON /SPACE/.  This  storage  is  used  only  during 
the  execution  of  the  subroutines  and  is  free  at  other  times  for  use  by 
other  routines. 
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4.2  SUMMATION  METHOD 

Another  method  of  calculating  G  is  to  sum  the  flux  from  different 
reflections  as  was  done  in  the  derivation  of  Eq.  (29).  From  Eq.  (15), 
if  Wi  is  set  equal  to  one,  then  Gjj  is  equal  to  Qy.  But  if  Wj  is  set  equal 

to  zero  for  j  not  equal  to  i,  then  by  definition  Qjj  is  equal  to  Qj  since 
the  flux  incident  on  the  jth  piece  had  to  come  from  the  ith  piece.  Thus, 
for  a  given  i,  if 

=  8{j  (56) 

then  letting  ^G  be  the  ith  row  of  G, 

iG  =  0'  (57) 

The  vector  Q  can  be  calculated  from  Eq.  (28).  The  first  term  of 
Eq.  (28)  is  calculated  from  Eq.  (26)  where  W  is  given  by  Eq.  (56). 
Succeeding  terms  are  obtained  from  Eq.  (27).  Letting  i  vary  from  one 
to  n,  G  can  be  calculated  a  row  at  a  time.  A  convergence  criterion  is 
easily  established  by  noticing  that  the  sum  of  the  elements  of  gives 
the  flux  remaining  after  k  reflections.  Since  a  unit  quantity  of  flux  was 
started,  the  summation  should  be  continued  until  the  quantity  of  flux 
left  is  small  (to  the  desired  accuracy)  when  compared  to  one. 

Subroutine  GSUM  was  written  to  calculate  G,  employing  the  method 
described  above.  This  routine  is  called  by  the  statement 

CALL  GSUM  (N,ND1M,R,F,G,NT,T0L) 

where  N,  NDIM,  R,  F,  and  G  are  the  same  as  for  GINV.  On  call,  NT 
contains  a  limit  to  the  number  of  terms  to  be  used  in  the  series.  On 
return,  NT  contains  the  number  of  terms  actually  used.  Accuracy  is 
controlled  by  the  argument  TOL.  Summation  is  continued  until  the  re¬ 
maining  flux  is  less  than  the  value  contained  in  TOL.  As  an  example, 
since  unit  flux  is  started;  given  a  value  of  0.  001  in  TOL,  summation 
continues  until  G  is  correct  to  three  significant  figures.  The  listing  of 
GSUM  is  given  in  Appendix  III. 

4.3  ITERATION  METHOD 

Multiplying  both  sides  of  Eq.  (22)  on  the  right  by  (I  -  RF)  and  trans¬ 
posing  the  negative  term  one  obtains 

G  =  F  +  GRF  (58) 

from  which  one  can  use  the  method  of  successive  approximations  to 
calculate  G.  That  is,  let 
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Gk  = 


F  +  G 


k  -  1 


RF 


for  k  =  0 
for  k  >0 


then 


G  =  lim  Gk 


(59) 


(60) 


Calculating  G  in  this  manner  would  require  one  two-dimensional  array 
to  store  Gk-1  and  another  two-dimensional  array  to  store  G^  and  thus 
would  use  more  storage  than  the  previous  methods.  However,  storage 
can  be  saved  by  calculating  G  one  row  at  a  time.  Writing  Eq.  (58)  in 
terms  of  the  elements,  it  becomes 


cu  -  F‘i  *  tl  Gi*"‘  «Ftk 


(61) 


Letting  ^G  and  ^F  be  the  ith  rows  of  G  and  F,  respectively,  Eq.  (61) 
becomes 

iGk  =  iF  -  iGk-1RF  (62) 


Equation  (62)  can  be  used  to  calculate  G  one  row  at  a  time. 


Using  Eq.  (58)  to  derive  an  expression  for  G,  one  obtains 

G  =  F  +  FRF  +  FRFRF  +  .  .  .  (63) 

When  calculating  G  a  row  at  a  time,  one  observes  that  each  term  of 
Eq.  (63)  is  the  contribution  attributable  to  the  corresponding  reflection. 
Thus  the  iteration  and  summation  methods  are  equivalent  except  that 
the  iteration  method  calculates  a  partial  sum  at  each  iteration,  where¬ 
as  the  summation  method  calculates  a  term  each  time  and  adds  it  to  the 
total.  Since  the  elements  of  ^G^  are  monotonically  increasing  with  k  and 
have  upper  limits,  convergence  is  ensured.  Equation  (39)  can  be  used 
as  a  convergence  criterion  to  test  for  any  desired  accuracy. 


It  is  interesting  to  note  that  a  formula  for  calculating  G  by  columns 
can  be  obtained  in  an  analogous  manner  except  using  Eq.  (33)  instead 
of  Eq.  (22).  Letting  Gj  and  Fj  be  the  jth  colums  of  G  and  F,  respec¬ 
tively,  the  corresponding  formula  is 

Gjk  =  Fjk  +  FRGjk_1  (64) 

Numerically,  calculating  G  by  columns  should  have  the  same  attributes 
as  calculating  it  by  rows;  however,  it  lacks  the  convenient  convergence 
criterion,  Eq.  (39). 
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Subroutine  GROW  was  written  to  calculate  G  by  iteration  a  row  at  a 
time.  This  routine  is  called  by  the  statement 

CALL  GROW(N,NDIM,R,F,G,NT,TOL) 

where  N,  NDIM,  R,  F,  and  G  are  the  same  as  for  GINV.  On  call  NT 
contains  a  limit  to  the  number  of  iterations  to  be  taken.  On  return  it 
contains  the  number  of  iterations  actually  taken.  Accuracy  is'controlled 
by  the  argument  TOL.  Iteration  is  continued  until  Eq.  (39)  is  satisfied 
within  the  given  tolerance.  Thus,  given  a  value  of  0.  001  in  TOL,  itera¬ 
tion  will  continue  until  G  is  correct  to  three  significant  figures.  The 
listing  of  GROW  is  given  in  Appendix  IV. 


SECTION  V 

CALCULATION  OF  THE  COEFFICIENTS 


5.1  DATA  MANAGEMENT 

In  the  process  of  taking  the  data  of  the  mathematical  model  and 
calculating  the  coefficients  of  Eqs.  (11)  through  (13),  three  tasks  are 
amenable  to  computer  programming  for  a  general  system.  These  are 
calculation  of  G,  conversion  of  units,  and  the  evaluation  of  the  equa¬ 
tions.  The  first  of  these  has  already  been  discussed.  The  subroutines 
UNIT  and  COEF  were  written  to  perform  unit  conversions  and  to 
evaluate  the  equations.  Before  describing  these  subroutines,  however, 
the  data  management  will  be  viewed. 

Table  I  tabulates  the  algebraic  variables,  the  FORTRAN  variables, 
and  the  dimensional  formulas  of  the  data  of  the  mathematical  model. 

For  the  moment  consider  only  the  storage  of  data.  It  can  be  seen  that 
the  first  eight  parameters  of  the  table  are  stored  in  the  same  array  P. 
This  was  done  solely  for  the  purpose  of  decreasing  the  length  of  the 
argument  list  of  the  subroutines.  The  first  eight  parameters  of  Table  I 
are  called  local  parameters,  since  they  give  the  properties  of  individual 
pieces.  Once  G  has  been  calculated,  the  form  factors  are  no  longer 
needed.  Thus,  if  this  calculation  is  done  before  the  conducting  areas 
and  distances  are  stored,  then  the  FORTRAN  array  F  can  be  equiva- 
lenced  to  the  variable  B.  Since  ay  and  dy  are  symmetric,  they  both 
can  be  stored  in  the  same  array  by  storing  ay  above  the  diagonal  and 
dy  below  the  diagonal,  replacing  the  form  factors.  Since  the  products 
diVici  occur  in  all  three  of  the  Eqs.  (11)  through  (13),  the  subroutine 
COEF  calculates  these  products,  replacing  the  ci.  The  coefficients,  Aj, 
replace  the  Pj8,  and  the  coefficients  By  replace  the  ay  and  dy.  Also, 
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the  coefficients  Cy  replace  the  Gy;  therefore  the  FORTRAN  array  G 
should  be  equivalenced  to  the  array  C.  With  this  arrangement  of  data, 
considerable  storage  economy  is  realized. 


TABLE  I 

DATA  REPRESENTATION  AND  DIMENSIONAL  FORMULAS 


Algebraic 

FORTRAN 

Dimensional  Formula 

di 

P(I,  1) 

M^i3 

Vi 

P(I,  2) 

T  3 

c.i 

P(I,  3) 

E4Mg1  ©e1 

*i 

P(I,  4) 

E7T81L92(0loLii)‘1 

«i 

P(I,  5) 

— 

ei 

P(I,  6) 

— 

Si 

P(I,  7) 

I  2 
^12 

Pf 

P(I,  8) 

El3T  [I 

aij 

B(I,  J),  J  >  1 

T  2 
l15 

dij 

B(I,  J),  J  <  I 

l16 

a 

SIG 

E17L18  ®  19 

diViCi 

P(I,  3) 

E  0  “ 1 

Fij 

F(I,  J)  S  B{I,  J) 

— 

Gy 

.  gg,  j)  *  ca,  j) 

— 

Ai 

P(I>  8) 

0T'1 

B(I,  J) 

T-1 

cij 

ca,  j) 

0  -3t-! 

5.2  UNIT  CONVERSIONS 

Subroutine  UNIT  was  written  to  perform  the  unit  conversions.  The 
subroutine  is  called  by  the  statement 


CALL  UNIT  (N,NDIM,SIG,P,B,FROM,TO) 


AEDC-TR-69-1  21 


where  briefly: 


N 

Contains  n 

NDIM 

Contains  the  dimension  of  the 
arrays  P  and  B 

SIG 

Returns  a  in  desired  units 

P 

Array  containing  local  parameters 
(See  Table  I) 

B 

Two-dimensional  array  containing 
conducting  areas  above  the  diagonal 
and  the  conducting  distances  below 

FROM 

Sixteen  element  integer  array 
indicating  units  of  input  data 

TO 

Nineteen  element  integer  array 
indicating  units  of  returned  data 

If  the  main  program  contains  the  statement 

DIMENSION  P  (100,8),  B(100,100) 

then  in  the  call  statement  above,  NDIM  would  contain  100.  The  dimen¬ 
sion  of  the  arrays  can  be  greater  than  or  equal  to  (not  less  than)  n.  The 
value  of  a  in  the  International  System  of  Units  (SI)  (Ref.  2), 

a  =  5.6697  x  10“*  — 

2  Oir* 

s  m  K 

is  stored  in  the  subroutine.  The  value  of  cr  as  indicated  in  the  array  TO 
is  returned  in  the  FORTRAN  variable  SIG.  When  UNIT  is  called,  F  and 
B  contain  data  in  the  units  indicated  in  FROM.  The  subroutine  replaces 
these  data  with  data  in  units  indicated  in  the  array  TO. 

The  dimensional  formulas  are  given  in  Table  I,  where  the  dimen¬ 
sions  are  indicated  as  follows: 


Symbol 

Dimension 

L 

Length 

M 

Mass 

T 

Time 

© 

T  emperature 

E 

Energy 

In  Table  I  the  input  data  have  subscripts  on  the  dimensional  symbols. 
These  subscripts  give  the  position  of  the  indicators  for  that  dimension 
in  the  arrays  FROM  and  TO,  Thus  the  fourth  elements  of  FROM  and  TO 
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contain  the  indicators  of  the  units  of  the  energy  dimension  of  the  specific 
heat  capacity.  It  is  noticed  that  there  are  two  length  dimensions  in  the 
dimensional  formula  for  the  thermal  conductivity.  The  first  derives 
from  the  conducting  area  and  the  second  from  the  temperature  gradient. 
The  formula  was  left  in  this  form  rather  than  simplified  since  it  is 
given  this  way  in  some  tables.  The  units  for  a  are  not  included  in  the 
array  FROM  since  its  value  is  not  input  to  the  subroutine,  but  rather 
is  stored  internally. 

The  indicators  of  the  various  units  and  the  conversion  factors  (ob¬ 
tained  from  Ref.  2)  are  tabulated  in  Table  II.  For  example,  values  of  4 
and  2  in  the  fourth  elements  of  FROM  and  TO,  respectively,  indicate 
that  the  unit  of  the  energy  dimension  of  the  specific  heat  capacity  data 
input  to  UNIT  is  in  Btu  and  it  is  desired  to  convert  it  to  ergs.  It  is  to 
be  noticed  that  a  one  in  each  element  of  TO  converts  all  units  to  SI.  The 
FORTRAN  listing  of  UNIT  is  given  in  Appendix  V. 


TABLE  II 

UNIT  INDICATORS  AND  CONVERSION  FACTORS 


Dimension 

Indicator 

Unit 

Conversion 

Factor 

Length 

1 

m 

1  m/m 

2 

cm 

0.  01  m/cm 

3 

in. 

0.  0254  m/in. 

4 

ft 

0.  3048  m/ft 

Mass 

1 

kg 

1  kg/kg 

2 

g 

0.001  kg/g 

3 

lb 

0.  4535933  kg/lb 

4 

slug 

14. 5939  kg/slug 

Time 

1 

s 

1  s/s 

2 

min 

60  s/min 

3 

hr 

3600  s/hr 

Temperature 

1 

°K 

1  °K/°K 

2 

°R 

0.  5555556  °K/°R 

Energy 

1 

J 

1  J/J 

2 

erg 

1  x  10"?  J /erg 

3 

cal 

4.  184  J/cal 

4 

Btu 

1054.  35  J  /Btu 

Note:  s  =  seconds 
J  =  joules 
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5.3  EVALUATION  OF  EQUATIONS 

Subroutine  COEF  was  written  to  evaluate  Eqs.  (11)  through  (13). 
It  is  assumed  that  all  data  are  in  consistent  units.  The  subroutine  is 
called  by  the  statement 


CALL  COEF  (N,NDIM,SIG,P,B>C) 

where: 

N  Contains  n 

NDIM  Contains  dimension  of  P,  B,  and  C 

SIG  Inputs  the  value  of  a  in  units  consistent 

with  the  rest  of  the  data 

P  Array  which  inputs  the  local  parameters, 

returns  the  product  dj  Vj  cj  in  place  of 
the  cj,  and  returns  the  coefficients  Aj 
in  place  of  the  Pi8 

B  Two-dimensional  array  which  inputs  the 

conducting  areas  above  the  diagonal  and 
the  conducting  distances  below  the  diag¬ 
onal;  returns  the  B  coefficients 

C  Two-dimensional  array  which  inputs  G  and 

returns  the  C  coefficients 


The  argument,  NDIM,  is  the  same  as  for  UNIT  except  that  it  applies 
also  to  the  array  C.  When  evaluating  Eq.  (12),  it  is  assumed  that  pieces 
are  in  contact  if  and  only  if  their  conducting  areas  are  positive.  If  cj  is 
zero,  then  Ai,  the  ith  row  of  B,  and  the  ith  row  of  C  are  set  to  zero  (see 
Section  2.  2).  The  evaluation  of  the  equations  is  straightforward,  and 
the  argument  descriptions  above  are  self-explanatory.  A  FORTRAN 
listing  of  COEF  is  given  in  Appendix  VI. 


SECTION  Yl 

SOLUTION  OF  THE  DIFFERENTIAL  EQUATIONS 


6.1  NUMERICAL  METHOD 

The  numerical  method  used  to  obtain  solutions  of  Eq.  (10)  is  the 
Taylor  expansion  of  order  four  (Ref.  1).  Leaving  off  the  summation 
signs  and  letting  it  be  understood  that  there  is  a  summation  over  j  from 
one  to  n,  Eq.  (10)  and  its  derivatives  are: 
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T;  -  A5  +  B £ j  Tj  +  Cij  Tj4  (65 

Ti  =  Tj  (Bij  +  4  Cij  Tj1)  (66 

Tj  =  Tj  (Bij  +  4  Cij  Tj1)  +  12  Cij  Tj2  Tj2  (67 

Ti  =  Tj  (Bij  +  4  Cij  Tj3)  +  24  Cij  Tj  Tj  (Tj2  +  1.5  Tj  Tj)  <68 

By  the  Taylor  formula  with  a  remainder  (Ref.  3)  the  temperature 
can  be  written 

Ti(t)  =  Ti(tl)  +  Ti  (O  (t  -  g  +  -L  Ti  (tt)  (t  -  g» 


where 


+  Ti(o  (t  -  o*  +  ±  t ;  a*)  (t  -  t,r 


tt  <  t*  <  t 


The  errors  obtained  by  approximating  Eq.  (69)  by 

TjW  =  T i  (t»)  +  Ti  (tj  (t  -O  +  y  Ti  (O  (t  -  t,)2 

1  1 

+  7  Tj  (O  (t  -  t,)1  +  -  Tf  (tl)  (t  -  tt)4 
6  24 


Ei  =  -  Ti  (t*)l  (t  -  t,)4  (72) 

If  the  temperature  error  per  unit  time  is  to  be  less  than  Ac,  then  for 
the  time  interval  from  tj  to  t,  it  is  necessary  to  ensure  the  inequality 

Ei  <  Ac  ft  -O  (73) 

for  all  i.  It  is  known  that  the  derivative  of  the  temperatures  with 
respect  to  time  of  the  physical  system  described  by  Eq.  (10)  approach 
constants.  This  indicates  that  the  absolute  value  of  the  derivatives  of 
higher  order  decreases  monotonically  with  t  to  zero.  Therefore,  from 
Eq.  (72) 

Ei  <  |  Tj  (tt) !  (t  -  t,)4  {74) 


Therefore,  to  guarantee  Eq.  (74),  it  is  sufficient  to  ensure  for  all  i 

1  .  •••• 

9tI  Tj  (tj)  I  (t  -  t,)4  <  Ac  (t  -  q)  / 


(t  -  t,)  < 


/24ACV/1 

\PFM7 
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To  find  the  maximum  time,  tm,  for  which  the  Taylor  expansions  for  all 
n  temperatures  are  valid  one  must  use  the  maximum  value  over  i  for 
the  denominator  of  Eq.  (76),  so 


/j.4Ac_V 


Given  the  initial  temperatures  and  time,  the  Taylor  expansion  and 
tm  are  calculated.  If  the  time,  t,  at  which  the  temperatures  are  de¬ 
sired  is  less  than  tm,  then  the  Taylor  expansions  are  evaluated  at  t. 
Otherwise  the  Taylor  expansions  are  evaluated  at  tm  and  then  tj  is  set 
equal  to  tm  and  new  expansions  and  a  new  tm  are  calculated.  This  is 
repeated  until  t  is  less  than  tm.  The  flow  chart  of  this  method  is 
shown  in  Fig.  3. 


Fig.  3  Flow  Chart  of  Numericol  Method  for  Obtaining  the  Solutions  to  Eq.  (10) 
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6.2  SUBROUTINE  SOLU 

Subroutine  SOLU  was  written  employing  the  numerical  method  de¬ 
scribed.  The  subroutine  is  called  by  the  statement 

CALL  SOLU  (N,NDIM,AC,Tl,T,TM,TEMP,Q,A,B,C) 

where 

N  Contains  n 

NDIM  Contains  the  dimension  of  Q,  A,  B,  and  C 

AC  Contains  the  allowable  temperature  error 

per  unit  time 

T1  On  first  call  contains  the  initial  time:  returns 

with  the  base  time  for  the  Taylor  expansions 
and  must  be  left  unchanged  for  subsequent 
calls. 

T  Contains  the  time  at  which  the  temperatures 

are  desired, 

TM  On  the  first  call  contains  the  initial  time; 

returns  with  tm  and  must  be  left  unchanged 
for  subsequent  calls 

TEMP  Returns  the  temperatures  at  the  time  con¬ 
tained  in  T 

Q  On  first  call  Q(l,l)  must  contain  the  initial 

temperature;  on  return  Q  contains  the  Taylor 
expansions  of  the  temperatures  and  must  be 
left  unchanged  for  subsequent  calls 

A  One-dimensional  array  containing  the  A 

coefficients 

B  and  C  Two-dimensional  arrays  containing  the 
B  and  C  coefficients,  respectively 

In  the  main  program,  if  the  arrays  are  dimensioned 

DIMENSION  TEMP  (70),Q(5,70),A(70),B(70,70),C(70,70) 

then  in  the  call  statement  above,  NDIM  would  contain  70.  The  dimen¬ 
sion  of  the  arrays  can  be  greater  than  or  equal  (not  less  than)  n.  Note 
that  if  the  local  parameters  are  no  longer  needed,  one  can  save  storage 
by  the  statement  (assuming  that  NDIM  is  the  same  for  COEF  and  SOLU) 

EQUIVALENCE  (P(1),Q(1)),  (P(1,7),TEMP  (1)),  (P(1,8),A(1)) 

The  FORTRAN  listing  of  SOLU  is  given  in  Appendix  VII. 
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SECTION  VII 
EXAMPLE  PROBLEMS 


* 


The  following  example  problems  are  to  serve  as  a  guide  for  users 
of  the  computer  programs.  The  problems  were  chosen  such  that  solu¬ 
tions  were  available  so  the  results  could  be  checked;  thus  they  serve  as 
check  problems  for  the  programs.  Although  the  examples  are  simple, 
the  real  usefulness  of  the  programs  is  for  much  more  complicated 
problems. 


7.1  ISOTHERMAL  RADIATING  OBJECT  WITH  A  HEAT  SOURCE 


The  differential  equation  for  an  isothermal  radiating  object  with  a 
heat  source  is 

dVcT  =  P"  -  SfffT4  (78) 

which  can  be  put  in  the  form 

T  =  A  +  CT4  (79) 


It  is  seen  that  the  steady-state  temperature,  T,,,,  is 


and  the  initial  time  derivative  of  the  temperature  is 

T0  =  A  +  CT04 

In  these  terms  the  solution  is 


(80) 

(81) 


A  program  was  written  using  SOLU  to  solve  Eq.  (79)  and  compare  the 
solution  with  Eq.  (82).  Some  of  these  solutions  are  shown  in  Fig.  3. 
The  accuracy  was  better  than  that  specified  by  the  FORTRAN  variable, 
AC,  and  the  error  was  much  too  small  to  be  seen  in  Fig.  4.  The 
FORTRAN  listing  of  this  program  is  given  in  Appendix  VIII. 
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7,2  CONDUCTING  RADIATING  SPHERICAL  SHELL 

Consider  a  spherical  shell  of  thickness  b  and  radius  a  shown  in 
Fig.  5.  The  emissivity  of  the  inside  surface  is  e  *  and  the  outside  sur¬ 
face  is  eO,  The  sphere  is  exposed  to  parallel  radiation  of  uniform  cross 
section  of  intensity  qs  for  which  the  absorptivity  of  the  surface  is  as. 
There  is  vacuum  inside  and  outside  the  sphere  and  it  is  assumed  that 
b  «  a.  The  heat -transfer  equations  for  this  system  have  been  derived 
in  Ref.  4,  and  a  numerical  method  of  solution  for  the  steady-state 
equations  is  given.  The  numerical  method  described  in  this  report  was 
applied  to  this  system  to  obtain  the  steady-state  solution  as  an  illustra¬ 
tion  of  the  application  of  the  method. 


Fig.  5  Conducting  Radiating  Spherical  Shell 
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Dividing  the  system  into  n  pieces  one  obtains  the  following  equa¬ 
tions: 


X  =  COB  6 

Ax  =  =  -sin  6  A  6 

n 


A0 

dij 


a  A  x 


2rrab  sin  6 


Pj®  =  (asqs)  (2ffosin  6)  (-  Aasintf) 
=  2  rra1  aa  qs  x  A  x 
Vj  =  2  n  aJ  b  Ax 


It  will  be  assumed  that  a 1  =  e*  and  a®  =  eO,  but  note  that  the  interior 
radiating  surface  of  a  piece  is  different  than  the  exterior  surface,  which 
is  contrary  to  the  assumptions  of  the  mathematical  method.  One  could 
divide  the  sphere  into  two  spherical  shells,  but  this  would  require  twice 
as  many  pieces.  This  can  be  avoided  by  defining  effective  emissivities 
and  form-surface  factors  in  the  following  manner.  The  radiating  area 
will  be  the  sum  of  the  internal  and  external  areas,  so 

S;  =  4  n  aJ  A  x 


The  effective  emissivity,  €  ',  will  be  defined  so  that  the  rate  radiant 
energy  is  emitted  remains  unchanged,  so 


Sif'alV  =  —  f^T.4  +  —  (°aT|4 
2  l  2  1 


or 


The  effective  form-surface  factor,  G^j,  will  be  defined  so  the  heat 
transfer  by  radiation  between  two  pieces  remains  unchanged,  so 


or 


a' Sj  e'  a  TV  G^j  =  a1  Si  el  a  Ty  Gjj 


Gij 


The  effective  form-surface  factor  is  thus 


29 


A  EDC-T  R-69- 121 


In  Ref.  4  it  is  proved  that  the  steady-state  temperature  distribution 
depends  only  on  the  two  parameters 


So  that  solutions  of  the  two  different  methods  could  be  compared,  the 
following  numerical  values  were  assigned  to  the  system  data  to  obtain 
the  specified  Nr  and  E: 

a  =  l 
£°  =  1 

aeqB  =  4(7 
b  =  a 


£>  =  E 


A  program  was  written  using  subroutines  COEF  and  SOLU  to  ob¬ 
tain  the  solutions.  A  FORTRAN  listing  of  the  program  is  given  in 
Appendix  IX.  Three  solutions  are  shown  in  Fig.  6.  Letting  n  =  20,  the 
solutions  by  the  two  different  methods  are  the  same  to  better  than  three 
significant  figures; 

7.3  CLAUSING  FACTOR  CALCULATIONS 

The  formula  for  the  form-surface  factor  was  based  on  the  assump¬ 
tion  of  diffuse  reflection.  Thus  the  use  of  the  formula  is  not  confined  to 
radiation,  but  can  be  used  to  calculate  flux  transfer  of  any  type  provided 
the  diffuse  reflection  assumption  is  valid.  In  fact,  the  method  has  been 
used  for  molecular  flux  as  reported  in  Ref.  5.  To  illustrate  the  method, 
subroutine  GINV  was  used  to  calculate  Clausing  factors.  The  Clausing 
factor  is  used  in  vacuum  technology  to  calculate  the  molecular  conduct¬ 
ance  of  tubes  and  is  defined  as  the  ratio  of  the  rate  at  which  gas  leaves 
the  outlet  of  a  tube  to  the  rate  at  which  it  is  incident  on  the  inlet.  To 
employ  the  method  the  tube  was  divided  into  n-2  equal  pieces  and  the 
inlet  and  outlet  were  each  counted  as  a  piece.  The  form  factors  for  the 
system  for  various  ratios  of  length  to  the  radius  of  the  tube,  L/r,  were 
calculated.  Each  piece  of  the  tube  was  assigned  a  reflectivity  of  one, 
and  the  inlet  and  outlet  were  assigned  a  reflectivity  of  zero.  It  can  be 
seen  that  the  Clausing  factor  is  just  Gin.  The  results  obtained  were 
compared  with  a  table  given  in  Ref.  6  and  are  shown  in  Fig.  7.  The 
FORTRAN  listing  of  the  program  used  is  given  in  Appendix  X. 
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APPENDIX  I 

FORTRAN  LISTING  OF  GINV 


SUBROUTINE  GINVI N, NDI M,R , F, G ,DET > 

C  CALCULATION  OF  G  BY  MATRIX  INVERSION  D.C.  TODD  10-21-68 
C  USES  4N  BYTES  OF  COMMON  /SPACE/.  RESERVES  400  BYTES. 

REAL  R(NDIM),FlNOIM,NDIM),G  { NDI M, NDI M ) , DET  ,V(100) 
COMMON  /SPACE/  V 
DO  1  1=1, N 
DO  1  J=1 ,N 

1  G( I , J )=-R ( I >*F (  I,J) 

DO  2  1=1, N 

2  G(I,I)=G(I,I>+1. 

CALL  INVERT(N»NDIM,G»DET> 

IFIDET.EO.O.) RETURN 
DO  6  J=1 ♦ N 
DO  3  1  =  1, N 

3  V ( I )=G< I,J> 

DO  5  1  =  1, N 
SUM=0. 

DO  4  K=1,N 

4  SUM=SUM+F< I ,K)*V(K> 

5  G( I  * J )  =  SUM 

6  CONTINUE 
RETURN 
END 
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APPENDIX  II 

FORTRAN  LISTING  OF  INVERT 


SUBROUTINE  INVERT  IN , IOI M.A.DET ) 

C  MATRIX  INVERSION  O.C.  TOOO  11-9-67 

C  USES  4N  BYTES  OF  COMMON  /SPACE/.  RESERVES  400  BYTES. 
INTEGER‘2  IRCI2.100) 

COMMON  /SPACE/  IRC 
DIMENSION  A ( IOIMv IOIMI 
C 

DET=1. 

00  9  K=L,N 
C  FIND  PIVOT 
ABSP-O. 

DO  L  I=K,N 
00  1  J-K .  N 
ABSA=  A9S I A  < I • J  > ) 

IFIA8SA.LE.ABSPIG0  TO  1 

IRC  1 1 ,X )  =  I 

1RCI2,X1=J 

ABSP=ABSA 

PIVOT-AI I,J) 

1  Continue 

C  OET  IS  PROOUCT  OF  PIVOTS.  IF  PIVOT  IS  ZERO.  THEN  A  I 
DET-OET-PIVOT 
I F ( PIVOT • EO.O. )GD  TO  15 
C  INTERCHANGE  ROWS  X  AND  I 
I  - 1 RC  ( 1 .  K  ) 

IF  1 1 «E0 .K )GD  TO  3 

DET=-OET 

DD  2  J  =  1.N 

SAVE-AIX.JI 

AII.JIMII.JI 

2  A(I,J)=SAVE 

C  INTERCHANGE  COLUMNS  X  AND  J 

3  J=IRC(2,X> 

IFIJ.EO.OGO  TO  5 
DET— DET 

00  4  I-l.N 
SAVE"A ( I ,K) 

All  .Kl-AI I.J1 

4  A(  I, Jl-SAVE 

C  OIVIDE  COLUMN  X  BY  MINUS  PIVOT 

5  R-l. /PIVOT 
DO  6  I-l.N 

6  All, XI—  R*At  I, XI 
C  REDUCE  MATRIX 

DO  20  I-l.N 

I F I I . E0.K1G0  TO  20 
00  7  J-l.N 

IFIJ. EQ.X1G0  TO  7 

A( 1,J)=A(I,J)+A(I,X>»A<X,JJ 

7  CONTINUE 
20  CONTINUE 

C  OIVIDE  ROM  X  BY  PIVOT.  AIK.KI-R 


DO  8  J  =  1,N 
B  A(K,J)-R*A(K,J) 
AIX.XI-R 
9  CONTINUE 
C  REARRANGE  MATRIX 
X=N*1 

10  K-X-l 

J= IRC  1 1 .X  I ' 
IFIJ.E0.K1G0  TO  12 
DO  11  1=1, N 
SAVE-A ( 1 ,K I 
A( I.KI-AI I,J) 

11  A( I, Jl-SAVE 

12  1= IRC (2 ,X I 

I F I I.EO.KIGO  TO  14 
DO  13  J=1 »N 
SAVE-A IX , J I 
A<X,X1=A(I,J! 

13  All, Jl-SAVE 

14  IFIK.NE.DGO  TO  10 
C 

15  CONTINUE 
RETURN 
ENO 
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AEDC.TR.49.121 


SUBROUTINE  GSUM ( N t NOI M,R ,F ,G ,NT ,TOL ) 

C  CALCULATION  OF  G  BY  SUMMATION  O.C.  TODD  10-21-68 
C  USES  8N  BYTES  OF  COMMON  /SPACE/.  RESERVES  800  BYTES. 
REAL  R{NDIM),F(NDIMtNDIM) tG (NOI Mt NDIM) ♦ V( 100t 2 ) 
COMMON  /SPACE/  V 
MM=0 

00  15  1=1 tN 
DO  10  J=1,N 
G( I  »  J )=F( I ,J> 

10  VIJ,1)=RIJ)*F(I,J) 

Ll  =  l 

L2=2 

DO  13  M=1,NT 
ERR=0. 

00  12  J  =  L  tN 
SUM=0. 

DO  11  K=1,N 

11  SUM=SUM+V(K,L1 )*FIK,J) 

G(  I  » J )=G I  I ,J)+SUM 

V ( J  tL2 )  =  R{J)*SUM 

12  ERR=ERR+V  C  J  »L2 ) 

I F ( ERR.LE.TOL  >G0  TO  14 

L3=L2 

L2  =  L1 

13  L1=L3 
M=NT 

14  IFIMM.LT.M)MM=M 

15  CONTINUE 
NT=MM 
RETURN 
END 
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FORTRAN  LISTING  OF  GROW 


C  CALCULATION  OF  G  BY  ROW  ITERATION  D.C.  TODO  10-21-TS8 
C  USES  8N  BYTES  OF  COMMON  /SPACE/.  RESERVES  800  BYTES. 
SUBROUTINE  GROW ( N, NDI M , R , F ,G ,NT » TOL ) 

REAL  R(NDIM) ,F(NDIM,NDIM) ,G(NDIM,NDIM),V( 100,2 ) 

COMMON  /SPACE/  V 

MM=0 

DO  6  1=1 ,N 
00  1  J=1,N 

1  VI Jtl)=F( IVJ) 

L  1  =  1 

L2=2 

00  4  M=1 , NT 
ERR=0 • 

DO  3  J  =  1,N 
SUM=0 . 

DO  2  K=1,N 

2  SUM=SUM+V(K,L1)*R(K)*F<K,J> 

V(J,L2)=F(I,J)+SUM 

3  £RR=ERR+V(J*L2)*(1.-R(J)) 

L3=L2 

L2=L  1 
L  1=L3 

IFI1.-ERR.LE.T0DG0  TO  5 

4  CONTINUE 
M=NT 

5  IFIMM.LT .M)MM=M 
DO  6  J=1,N 

6  G(I,J)=V(J,L1) 

NT=MM 

RETURN 

END 
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APPENDIX  V 

FORTRAN  LISTING  OF  UNIT 


SUBROUTINE  UNIT ( N,NDIM , S IG ,P ,B, FROM,  TO ) 

C  CONVERT  UNITS  OF  SYSTEM  DATA  D.C.  TODD  10-24-68 
REAL  PINDIM,8),B(N01M,NDIM), 

L  L(4)/l.,. 01, .0254, .3048/, 

M  M (4 )/l.,. 001,. 4535933, 14. 5939/, 

S  S(3)/l.,60.,3600./, 

T  T(2)/l.,. 5555556/, 

E  E ( 4 ) /I . , 1 . E-7, 4. 1 84, 1054. 35/ 

INTEGER  FROM! 16 ) ,T0( 19) 

S IG=5. 6697E-8 

F0= I  Ml  FROM (1>)/M(T0I1>) )*<L(T0(2> >/L( FROM! 2) ) )**3 
FV=(L(FROM(  3)  >/L(T0(3)>>**3 

FC=E  I FROM (4))*M(T0I5))*T(T0(6)>/(EIT0(4) ) *M( FROM! 5 ) )*T( FROMI 6 ) ) ) 
FK=E  I FROM  17) )*L  I FROM( 11) ) *S I  TO  1 8 ) )*T ( T0( 10 ) )* ( L ( TO ( 9 ) >/L( FROM! 9) ) ) 
K**2/(EIT0(7)  >*UTO<  11 )  >*S<  FROM! 8)  >*T(  FROMI  10)  )  ) 

FS= I L I FROM( 12 ) ) /L ( TO ( 1 2 ) ) ) **2 

FP=E ( FROM  1 13) )  *S  I T0( 14) )/ (E( T0(13) >  *  S I FROMI 14) >) 
FA=(L(FROM(15))/L(TO(15> ) >**2 
FN=L ( FROM  116) >/L(T0(16> ) 

FG=(E( FROMI 17) )/E(TO( 17) ) >*(L I  TO < 18 > > /L ( FROMI 18) ) )**2 
G*(T(TOI 19) )/T (FROMI 19) ) )**4 
00  1  1*1 ,N 
PI  1  , 1 ) =FD*P (1,1) 

P(I,2)*FV*P( 1,2) 

PI I,3)=FC*P( 1,3) 

PI  1 ,4)*FK*P( 1,4) 

PI I,7)«FS*P( 1,7) 

1  P I  1 , 8 )=FP*P (1,8) 

NM1=N— 1 

DO  2  1=1, NM1 
I  PI* I +1 
DO  2  J=IP1,N 
B  I  I , J )=FA*B 1 1 , J  ) 

2  B(J,I)*FN*B(J,  I) 

S1G=FG*SIG 

RETURN 

END 
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SUBROUTINE  COEF ( N .NO IM , S 1 G , P , B ,C ) 

C  COEFIC IENTS  FOR  HEAT  TRANSFER  DIFFERENTIAL  EQUATIONS 
C  O.C.  TODO  10-23-6B  ST3OO2-YO0 
C  P  INPUTS  D,V,C,K, ALPHA, EPSILON. S,  AND  P 
C  RETURNS  DVC  IN  C  AND  A  IN  P 

C  B  INPUTS  AREAS  ABOVE  THE  DIAGONAL  ANO  DISTANCES  BELOW 
C  RETURNS  WITH  8  COEFIC'lENTS 
C  C  INPUTS  G  AND  RETURNS  WITH  C  COEFFICIENTS 

REAL  P<NDIM,B),B(NDIM,NDIM),C(NDIH,NOIM) 

DO  2  1=1, N 

D=PCI,l‘)*PCI,2)*PCI,3) 

P( 1,31=0 

IFJD. NE.O.IGO  TO  1 
PI  1,81=0. 

GO  TO  2 

1  P(I,ai=P(I,81/D 

2  CONTINUE 
DO  7  1=1, N 
D=P( 1,3) 

IFIO. NE.O. 1G0  TO  24 

B ( I T I  1=0. 

C (  I, I  1=0. 

GO  TO  7 
24  SUM=0. 

DO  6  K=1,N 
IF(K-I14,6,3 

3  II*I 
JJ=K 

GO  TO  5 

4  1 1=K 
JJ  =  I 

5  IFIBIII, JJ1.EQ.0. 1G0  TO  6 
SUM=SUM+P<K,4)*BIII, J J  > / <  B  <  J J , 1 1 >  *  I  PC I,4)+PIR,4)) I 

6  CONTINUE 

B( I » I )=-2.*P  < I»41*SUM/D 

C Cl, I 1=SIG*P (I, 6 )*P< I, 7)4 IPI! ,5) *C (I, I >-l. )/D 

7  CONTINUE 
NM1=N-1 

00  13  1=1, NM1 
DI=P( I ,3) 

IP1=I+1 
DO  12  J=IP1,N 
DJ=P ( J , 3 ) 

IF(B(I,J).EQ.O.)B(J,I )=1. 

11  E=2.*PC I,4)*P( J,4)*B( I ,J)/(B( J,I >*(P(I,4)+P( J,4) )) 
IFCOI.NE.O. )G0  TO  20 
B ( I , J 1=0. 

CIJ=0. 

GO  TO  21 

20  B I  I , J 1 =E/OI 

CIJ=SIG4P( I ,5)*C( J, I )*P< J,7)*PI J,6)/DI 

21  IFCDJ.NE.O. ) GO  TO  22 


BU,  1 1=0. 

C( J, I  1=0. 

GO  TO  23 

22  B ( J , I  1 =E/0 J 

C(  J,I  )=SIG4P(J,5)*C(  I  ,J)*P(I  ,7)*PU,6)/DJ 

23  C ( I  ,  J 1 =C  I J 

12  CONTINUE 

13  CONTINUE 
RETURN 
END 
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SUBROUTINE  SOLU ( N, NOI M, AC , T1 ,T, TM, TEMP , P, A, B, C ) 

C  SOLUTION  TO  HEAT  TRANSFER  DIFFERENTIAL  EQUATION 
C  D.C.  TODD  10-25-68  ST8002-Y0O 

C  AT  FIRST  CALL  P ( 1  *  I >  INPUTS  INITIAL  TEMPERATURES  AND  TM=T1 
C  CONTENTS  OF  P,  Tl,  AND  TM  MUST  BE  LEFT  UNCHANGED  FOR  SUBSEO0ENT  CALLS 
REAL  TEMP(1),P(5,1)»A(1>,B(NDIM,1>,C(NDIM,1) 

1  IF(TM.GE.T)GO-TO  2 
X=TM 

K=0 

GO  TO  3 

2  X=T 
K= 1 

3  DX=X-T 1 

DO  4  1  =  1  ,N 

4  TEMPI  I )  =  P 1 1  *  I >+DX*( PI2, 1 >+DX*(P(3, I >+DX*(P(4, I )+DX*P(5, I ) ) ) ) 
IFIK.EQ.l > RETURN 

T  1  =  TM 
P5M=0 • 

DO  5  1=1, N 
PI  1 , 1 )=TEMP( I  ) 

P (2, 1 )=A  1 1 ) 

P 1 3 , 1 )=0 . 

P 1 4, 1 )=0. 

5  PI5,I)=0. 

DO  6  1=1, N 
DO  6  J=1,N 

6  P(2,I)=P(2,I)+BI I, J)*P(1,J)+CI I,J)*PI1, J)**4 
DO  7  1=1, N 

DO  7  J=1,N 

7  PI3»I )=P(3,I)+P(2,J)*(B(I,J) +4.*C( I,J)*PI1,JI**3) 

DO  8  1=1, N 

DO  8  J=1,N 

8  P(4,I  )  =  P(4,I  >  +  P{3,J>*(BU,J>+4.*C(I,  J)*P(1,  J)**3) 

1  +12. *C 1 1 , J »*P 1 1, J )**2*P 1 2 , J )**2 

DO  10  1  =  1, N 
DO  9  J=1,N 

9  P I  5 , I )  =  P( 5 , 1 >+PI4, J )* (B( I , J )+4. *C( I , J )*P 1 1, J )**3) 

1  +24. *C 1 1 , J )*P 1 1 , J )*P 1 2 , J )*IP 1 2, J )**2+l,5*P 1 1 , J ) *P 1 3, J ) ) 

AP5=ABSI PI5, I ) ) 

IFIAP5.GT  *P5M)P5M=AP5 

10  CONTINUE 
P5M=P5M/24. 

DO  11  1=1, N 

P 1 3 , 1 )  =  ,5*P 13,1) 

P 14, I )=. 1666 667* P (4,1) 

11  P I  5 , 1 )  =  .04166667*P 15,1) 
iFlPSM.GT. 1. E-66 )G0  TO  12 
TM=1.E16 

GO  TO  1 

12  TM=T1+(AC/P5M)**. 3333333 
GO  TO  1 

END 
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FORTRAN  LISTING  OF  PROGRAM  TO  COMPARE  SOLUTIONS  OF  RADIATING  OBJECT 


C  COMPARISON  OF  ANALYTICAL  ANO  NUMERICAL  SOLUTIONS  OF  ISOTHERMAL 
C  RADIATING  OBJECT  WITH  HEAT  SOURCE  D.C.  TODD  12-16-68  AEA00071 

REAL  P<  5) , TI*8,CONST*8 

1  READ(5,1001 , END=900 )N,TO,TI , TOD, AC 
IF (TI .NE.O. )GO  TO  2 

A=0. 

C=T0D/T0**4 

D=1./(3.*C) 

CONST =1 .DO/ T0*#3 
GO  TO  3 

2  A=TOD/!1.-<TO/TI)**4) 

C=-A/TI**4 

D=T I / (4.*A ) 

CONST =DLOG ( (TO+T I ) /DABS < TO-TI ) >+2.*DATAN( TO/TI ) 

3  S1=0. 

SM=0. 

P( 1 )=T0 

DT= (T I-TO )/ (N+l ) 

T=TO 

WRITE! 6,1000) 

WRITE (6,1005) 

WRITE! 6,1003 )T0,TI,T00,AC,A,C,0, CONST 
WRITE ( 6 ,1006 ) 

DO  6  1=1, N 

t=t+dt 

I F ! TI • NE.O. )G0  TO  4 
S=D* (CONST-1. 00/T**3) 

GO  TO  5 

4  S=D* ( DLOG ( ( T+T I ) /DABS! T— TI ))+2.*DATAN(T/TI )-CONST> 

5  CALL  S0LU(1,1,AC,S1,S,SM,TEMP,P,A,0.,C) 

ERR=T-TEMP 

6  WRITE(6,1003)S,T, TEMP, ERR 
GO  TO  1 

900  WRITE ( 6, 1004) 

STOP 

1000  FORMAT! 'ID. C.  TODO  AEA00071') 

1001  FORMAT! 12, E10 .0 ,5E12 .0 > 

1003  FORMAT (1P11E12.4) 

1004  FORMAT! *0THE  END‘/1H1> 

1005  FORMAT ( 1H04X2HT010X2HTI 10X3HT0D9X2HAC11X 1HA1 1X1HC1 IX 1HD9X5HC0NST / ) 

1006  FORMAT! 1H05X1HS 1 1XIHT9X4HTEMP9X3HERR/ ) 

END 
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FORTRAN  LISTING  OF  PROGRAM  TO  OBTAIN  SOLUTION 
TO  CONDUCTING  RADIATING  SPHERICAL  SHELL 


C  CONDUCTING  RAO  I  AT  I NG  SPHERICAL  SHELL  O.C.  TODO  I2-16-6B  AEA00072 

REAL  PI60.B),0!5,601,A(6D1 ,B!6D, 6D1 ,C( 60,60) .TEMPt 60) .NR, 

*  P 1/3. 1*1593/, S1G/5.6697E-B/.TURE I 6D1 
EQUIVALENCE  I P ( 1 ) ,0< 1 ) ) , IP  1 1 ,7) , TEMP  ( 1 1 ) , t PI  1 , B ) , A  Cl) ) 

1  RE AO (5,1001, EN0=90D 1NH,E  »NR 
NR  IT  E (6, 1000 1 
WRITE(6,1005)E,NR 

N=2*NH 
DX-2./N 
DS=*.*PI*DX 
DV»2.*PI=SIG*DX 
TK  =  1 ./NR 
6P=.5*I£+1.  ) 

AQ=B.*P1«SIG*DX 
x=-l.-.5*0X 
DD  2  1=1, N 
X=X*OX 
P(1,1>*1. 

P( l,2)=DV 
PII,31=TK 
P(I,*)=TK 
PI  1 ,5 1=EP 
P  1 1  •  6  )*EP 
PCI. 71-OS 

2  P ( I , 8 1=A0*X 
DO  3  1=1 ,NH 

3  P(I,B>=0. 

G=.25*E»0X/EP**2 
DD  *  1=1, N 

00  *  J=1 ,N 
B ( I , J 1=0. 

*  C ( I , J )=G 
x=-l. 

DO  5  J=2,N 

X-X+DX 

I-J-l 

ST=(1.-X»X1**.5 

BCI,JI»2.*PI*S1G»ST 

5  8 ( J , I 1=DX/5T 

CALL  C0EF(N,60,SIG«P,B,C1 
T1=0. 

T-l.E-10 

TH=0. 

00  6  1=1, N 
TURE ( I  1  =  1. 

6  0(1.1 1=1. 

WRITE ( 6, 1006  ) 

LAST=0 

I PRT=20 

DO  9  IT-1.1D0 

CALL  S0LUIN,60,.2,T1,T,TM,TEMP,0,A, B,C) 

DMAX=0 • 


DO  7  1=1, N 
OIFP=TENP(I)-TURE(I> 
IF(ABS(DIFF).GT.ABS(OMAX) )OMAX=DIFF 
7  TUREI I )=TEMP(I 1 

WRITEI6.1 002 )IT,T,OHAX 
IF(T .LT.I.E151G0  TO  B 
I PRT= I T 
LAST-1 

B  IFIIPRT.NE.ITIGO  TO  9 
I PRT= I PRT+20 

HRITE(6,1DD3)(TEMP( | ),I=1,N) 
IFILAST.NE.01G0  TO  10 
9  T=1.D01»TM 
ID  GO  TO  1 
90D  WRITE! 6,100*1 
STOP 

10D0  FORMAT 1< 10. C.  TODD  AEA00D71'! 

10D 1  F0RMAT(12,E1D.D,5E12.D> 

1002  FORMAT (I16.6X1P2E12.*) 

1003  FORMAT ( 1P1 1E12.A 1 
100*  FORMAT! >OTHE  ENO’/lHl! 

1005  FORMAT) >0E  =,1PE12.*,6X'NR  ='E12.*1 
IDO 6  FORMAT (1H  1 
ENO 
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FORTRAN  LISTING  OF  PROGRAM  TO  CALCULATE  CLAUSING  FACTORS 


C  CALCULATION  OF  CLAUSING  FACTORS  D.C.  TODD  1-3-69 
REAL  F ( 60 , 60 ) f A  <  2 ) , G ( 60 , 60 ) »  R I  60 ) 

WR IT6 ( 6  r 1000 ) 

1  READ(5tl001*  END* 900 )M*H 
CALL  TUBE (M»1.»H»A»F) 

N=M+2 

DO  2  I=2»N 

2  R( 11*1. 

R 1 1  )=0. 

R { N ) =0 • 

CALL  GINV(Nt60tR»F,G,DET) 

IFIDET.EO.O. )G<1,N)=2. 

WRITEI6,1003)M,H,G( 1,N> 

GO  TO  1 

900  WRITE ( 6 t 1004 ) 

STOP 

1000  FORMAT ( * 1AEA00075  D.C.  TODD'/) 

1001  FORMAT (I2yE22.0) 

1003  FORMAT ( I6»6X»1P2E12.4> 

1004  FORMAT ( • OTHE  END'/lHl) 

END 
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The  resulting  system  of  first-order  ordinary  differential  equations  is 
put  in  a  form  convenient  for  programming  a  numerical  solution  by 
defining  coefficients  in  terms  of  the  data  of  the  system.  Form-surface 
factors  are  a  generalization  of  form  factors  to  include  multiple  reflec¬ 
tions.  The  factor  is  defined  and  discussed.  A  formula  is  derived  for 
the  factor  in  terms  of  the  form  factors  and  the  reflectivities.  The 
form-surface  factor  relations  for  the  conservation  of  energy  and  the 
reciprocity  law  are  derived.  The  conditions  necessary  and  sufficient 
for  its  existence  are  derived.  Methods  of  calculation  are  discussed.  A 
FORTRAN  subroutine  is  presented  that  calculates  the  form-surface  factors 
from  the  form  factors  and  the  reflectivities. 
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